# ACCELERARE il codice controllando sempre se N<T o N>T e quindi fare X'X o XX'

#START  TIME VARYING  NUMBER OF FACTORS

library("mgcv")
library("parallel")

# Funzione per la generazione dei dati
{  
#  Dati:
#    - corr.cs: correlazione costante cross section
#    - corr.ts: parametro rho di un AR(1) per introdurre correlazione tra gli errori
#    - heterosk.ts : (qualunque valore) introdurre eteroshedasticita' nell'andamento temporale degli errori
  
  generazione_dati = function(N, T, r, theta, heterosk.ts = NA, corr.ts = NA, corr.cs = NA)
  {
    loadings = matrix(rnorm(N*r), nrow = N, ncol = r)
    factors = matrix(rnorm(T*r), nrow = r, ncol = T)
    c = loadings %*% factors
    e = matrix(rnorm(N*T), nrow = N, ncol = T)
    
    if (is.na(heterosk.ts) & is.na(corr.ts) & is.na(corr.cs))  
    {
      X = c + sqrt(theta) * e
      return(X)
    }
    
    if (!is.na(heterosk.ts))  
    {
      ts = diag(c(1, sqrt(2)), T)
      e = e %*% ts
      X = c + sqrt(theta) * e
      return(X)
    }
    
    if (!is.na(corr.ts) & is.na(corr.cs))  
    {
      e[,1] = e[,1]*sqrt(1/(1-corr.ts^2))
      for (t in 2:T) e[,t]= corr.ts*e[,t-1] + e[,t]
      X = c + sqrt(theta) * e
      return(X)
    }
    
    if (is.na(corr.ts) & !is.na(corr.cs))  
    {
      J = max(N/20, 10)
      e2 = e
      for (i in 1:N) e[i,] = (1-2*corr.cs)*e2[i,] +corr.cs*(sum(e2[(max(0, i-J)):i]) + sum(e2[i:(min(N, i+J))]))
#      cs = diag(0, N)
#      supp = c(rep(0, N-1-J), rep(corr.cs, J), 1, rep(corr.cs, J), rep(0, N-1-J))
#      for (i in 1:N) cs[i,] = supp[(length(supp) - i + 1 - N + 1):(length(supp) - i + 1)]
#      e = corr.cs %*% e
#     molto lento a lavorare con una matrice grande (80 secondi per generare 80 set di dati)
      X = c + sqrt(theta) * e
      return(X)
    }
    
    if (!is.na(corr.ts) & !is.na(corr.cs))  
    {
      e[,1] = e[,1]*sqrt(1/(1-corr.ts^2))
      
      J = max(N/20, 10)
      e2 = e
      for (i in 1:N) e2[i,] = (1-2*corr.cs)*e[i,] +corr.cs*(sum(e[(max(0, i-J)):i]) + sum(e[i:(min(N, i+J))]))
      
      for (t in 2:T) e[,t]= corr.ts*e[,t-1] + e2[,t]
      X = c + sqrt(theta) * e
      return(X)
    }
  }
}
# Funzione unica per tutti i criteri
{
  V = function(X, k)
  {
    N = nrow(X)
    T = ncol(X)
    lambda = sqrt(N) * as.matrix(eigen(X %*% t(X))$vectors[,1:k])
    factor = 1/N *(t(X) %*% lambda)
    out = 1/(N*T) * sum((X - lambda %*% t(factor))^2)
    return(out)
  }
  
  criteria = function(X, k, kmax, eps1 = NA, eps2 = NA)
  {
    N = nrow(X)
    T = ncol(X)
    out1 = V(X, k)
    out2 = V(X, kmax)
    
    out3 = k * ((N+T)/(N*T)) * log((N*T)/(N+T))
    PC_P1 = out1 + out2 * out3
    IC_P1 = log(out1) + out3
    
    out3 = k * ((N+T)/(N*T)) * log(min(N,T))
    PC_P2 = out1 + out2 * out3
    IC_P2 = log(out1) + out3
    
    out3 = k * log(min(N,T))/min(N,T)
    PC_P3 = out1 + out2 * out3
    IC_P3 = log(out1) + out3
    
    out3 = k * 2/T
    AIC1 = out1 + out2 * out3
    
    out3 = k * log(T)/T
    BIC1 = out1 + out2 * out3
    
    out3 = k * 2/N
    AIC2 = out1 + out2 * out3
    
    out3 = k * log(N)/N
    BIC2 = out1 + out2 * out3
    
    out3 = k * 2 * ((N + T - k)/(N*T))
    AIC3 = out1 + out2 * out3
    
    out3 = k * ((N + T - k)*log(N*T)/(N*T))
    BIC3 = out1 + out2 * out3
    
    out = c(PC_P1, PC_P2, PC_P3, IC_P1, IC_P2, IC_P3, AIC1, BIC1, AIC2, BIC2, AIC3, BIC3)
    
    if (is.na(eps1)) return(out)
    
      for (i in 1:length(eps2))
    {
      out3 = k * ( (log(N))^eps1 * N^eps2[i])/(T*sqrt(N)*sqrt(N))
      NEW1 = out1 + out3
      NEW2 = log(out1) + out3
      NEW3 = out1 + out2 * out3
      NEW4 = T/(T-k)*out1 + out3
      out = c(out, NEW1, NEW2, NEW3, NEW4)
    }
    return(c(out, out1))
  }
}
# funzione per costruzione tabella
{
  temp = function(count, N, T, r, theta, heterosk.ts, corr.ts, corr.cs, kmax, ncol, eps1, eps2)
  {
    X = generazione_dati(N, T, r, theta, heterosk.ts, corr.ts, corr.cs)
    temp2 = matrix(NA, nrow = kmax, ncol = ncol - 2)
    #kmax = min(kmax, T, N)
    for (k in 1:kmax) temp2[k,] = criteria(X, k, kmax, eps1, eps2)
    out = rep(NA, ncol - 2)
    for (l in 1:(ncol - 2)) out[l] = which(temp2[,l] == min(temp2[,l]))
    return(out)
  }
  
  
  tabella = function(r, theta, kmax, N.vec, T.vec, n.iter, eps1 = NA, eps2 = NA, heterosk.ts = NA, corr.ts = NA, corr.cs = NA)
  {
    t = Sys.time()
    PC.p1 = PC.p2 = PC.p3 = rep(NA, length(N.vec))
    IC.p1 = IC.p2 = IC.p3 = rep(NA, length(N.vec))
    AIC.1 = AIC.2 = AIC.3 = rep(NA, length(N.vec))
    BIC.1 = BIC.2 = BIC.3 = rep(NA, length(N.vec))
    table1 = data.frame(N.vec, T.vec, PC.p1, PC.p2, PC.p3, IC.p1, IC.p2, IC.p3, AIC.1, BIC.1, AIC.2, BIC.2, AIC.3, BIC.3)
    if (!is.na(eps1))
    {
      NEW.1.1 = NEW.2.1 = NEW.3.1 = NEW.4.1 = rep(NA, length(N.vec))
      NEW.1.25 = NEW.2.25 = NEW.3.25 = NEW.4.25 = rep(NA, length(N.vec))
      NEW.1.45 = NEW.2.45 = NEW.3.45 = NEW.4.45 = rep(NA, length(N.vec))
      NoPen = rep(NA, length(N.vec))
      table1 = data.frame(table1, NEW.1.1, NEW.2.1, NEW.3.1, NEW.4.1,  NEW.1.25, NEW.2.25, NEW.3.25, NEW.4.25, NEW.1.45, NEW.2.45, NEW.3.45, NEW.4.45, NoPen)
      rm(NEW.1.1, NEW.2.1, NEW.3.1, NEW.4.1,  NEW.1.25, NEW.2.25, NEW.3.25, NEW.4.25, NEW.1.45, NEW.2.45, NEW.3.45, NEW.4.45, NoPen)
    }
    rm(N.vec, T.vec, PC.p1, PC.p2, PC.p3, IC.p1, IC.p2, IC.p3, AIC.1, BIC.1, AIC.2, BIC.2, AIC.3, BIC.3)
    
    for(l in 1:nrow(table1))
    {
      print(l)
      temp1 = matrix(NA, nrow = n.iter, ncol = ncol(table1) - 2)
      temp3 = parLapply(cl, 1:n.iter, temp, N = table1[l, 1], T = table1[l, 2], r = r, theta = theta, heterosk.ts = heterosk.ts, corr.ts = corr.ts, corr.cs = corr.cs, kmax = min(kmax, table1[l, 1], table1[l, 2]), ncol = ncol(table1), eps1 = eps1, eps2 = eps2)
      for (j in 1:n.iter) temp1[j,] = temp3[[j]]
      table1[l, 3:ncol(table1)] = apply(temp1, 2, mean)
      #gc()
    }
    rm(l, j, temp1, temp3)
    print(Sys.time() - t)
    return(table1)
  }
  
    
 
  
}

no_cores <- detectCores() - 2
cl <- makeCluster(no_cores)
clusterExport(cl, c("generazione_dati", "criteria", "V"))

# temp3 = parLapply(cl, 1:10000000, sqrt)
# rm(cl)
gc()
memory.size(max = F)

# generate parameters from data - calibration

# import data using the 
Tmezzi=480/2
Tquarti=480/4
# T/4
{

delete_bal = apply(data[(Tmezzi+Tquarti+1):480,], 2, function(x){sum(is.na(x))})
delete_bal = which(delete_bal > 0)
N_bal=ncol(data[,-delete_bal])
X=data[(Tmezzi+Tquarti+1):480,-delete_bal]-(rf[(Tmezzi+Tquarti+1):480] - 1)%*%t(rep(1,N_bal)) # i : i+w-1
Xd=X-rep(1,Tquarti)%*%t(apply(X,2,mean))

meanX=apply(X,2,mean)


lambda_dgp =  sqrt(N_bal) * as.matrix(eigen(t(Xd) %*% (Xd))$vectors[,1:10])
factors_dgp= (1/N) *( Xd %*% lambda_dgp )

library(expm)

mean_lam = apply(lambda_dgp,2,mean)
sd_lam = sqrtm(cov(lambda_dgp))



mean_fac = apply(factors_dgp,2,mean)
sd_fac = sqrtm(cov(factors_dgp))


Ve=matrix(0,N_bal,10)
j=1
for (j in 1:10)
{ 
  e=Xd-factors_dgp[,1:j]%*%t(lambda_dgp[,1:j])
  Ve[,j]=diag(cov(e))
}

}

#T/2
delete_bal = apply(data[(Tmezzi+1):480,], 2, function(x){sum(is.na(x))})
delete_bal = which(delete_bal > 0)
N_bal=ncol(data[,-delete_bal])
X=data[(Tmezzi+1):480,-delete_bal]-(rf[(Tmezzi+1):480] - 1)%*%t(rep(1,N_bal)) # i : i+w-1
Xd=X-rep(1,Tmezzi)%*%t(apply(X,2,mean))

meanX=apply(X,2,mean)


lambda_dgp =  sqrt(N_bal) * as.matrix(eigen(t(Xd) %*% (Xd))$vectors[,1:10])
factors_dgp= (1/N) *( Xd %*% lambda_dgp )

library(expm)

mean_lam = apply(lambda_dgp,2,mean)
sd_lam = sqrtm(cov(lambda_dgp))



mean_fac = apply(factors_dgp,2,mean)
sd_fac = sqrtm(cov(factors_dgp))


Ve=matrix(0,N_bal,10)
j=1
for (j in 1:10)
{ 
e=Xd-factors_dgp[,1:j]%*%t(lambda_dgp[,1:j])
Ve[,j]=diag(cov(e))
}





# mio codice non parall
{
  M=3
  mm=1
  
  eps1 = 1e-8
  eps2 = 0.25*c(0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80) # 0.68 good for Tshalf
  
  set.seed(12)

  T=500
  N=N_bal#283
  Xsim=matrix(NA,T,N)
  Fsim1=rnorm(100,0,1)
  Fsim2=rnorm(200,0,1)
  Fsim2=matrix(Fsim2,100,2)
  Fsim4=rnorm(400,0,1)
  Fsim4=matrix(Fsim4,100,4)
  Fsim5=rnorm(500,0,1)
  Fsim5=matrix(Fsim5,100,5)
  Fsim8=rnorm(800,0,1)
  Fsim8=matrix(Fsim8,100,8)
  Fsim42=rnorm(400,0,1)
  Fsim42=matrix(Fsim42,100,4)
  Fsim22=rnorm(200,0,1)
  Fsim22=matrix(Fsim22,100,2)
  Fsim12=rnorm(100,0,1)
  
  Lsim1=rnorm(N,0,1)
  Lsim2=rnorm(N*2,0,1)
  Lsim2=matrix(Lsim2,N,2)
  Lsim4=rnorm(N*4,0,1)
  Lsim4=matrix(Lsim4,N,4)
  Lsim5=rnorm(N*5,0,1)
  Lsim5=matrix(Lsim5,N,5)
  Lsim8=rnorm(N*8,0,1)
  Lsim8=matrix(Lsim8,N,8)
  Lsim22=rnorm(N*2,0,1)
  Lsim22=matrix(Lsim22,N,2)
  Lsim42=rnorm(N*4,0,1)
  Lsim42=matrix(Lsim42,N,4)
  Lsim12=rnorm(N*1,0,1)
  
  kmax=10#10
  ncol=12+32+1 #12+4*eps2+1
  

 # eps2 = 0.65*c(0.25, 0.45) # 0.75 ok
  outt=matrix(NA,kmax,ncol)

  for (mm in 1:M)
  {
  # here simulate data of length T=480 and N=283
  
    # simulation
    for (t in 1:100)
      Xsim[t,] = t(meanX) + 
        (sd_fac[1,1]*Fsim1[t] + mean_fac[1]) %*% t(Lsim1*sd_lam[1,1] + mean_lam[1]) + 
        t(sqrt(diag(Ve[,1])) %*% rnorm(N,0,1))
    for (t in 101:200)
      Xsim[t,] = t(meanX) + 
        (Fsim4[(t-100),] %*% sd_fac[1:4,1:4] + t(mean_fac[1:4])) %*% t(Lsim4 %*% sd_lam[1:4,1:4] + rep(1,N) %*% t(mean_lam[1:4])) + 
        t(sqrt(diag(Ve[,4])) %*% rnorm(N,0,1))
    for (t in 201:300)
      Xsim[t,] = t(meanX) + 
        (Fsim8[(t-200),] %*% sd_fac[1:8,1:8] + t(mean_fac[1:8])) %*% t(Lsim8 %*% sd_lam[1:8,1:8] + rep(1,N) %*% t(mean_lam[1:8])) + 
        t(sqrt(diag(Ve[,8])) %*% rnorm(N,0,1))
    for (t in 301:400)
      Xsim[t,] = t(meanX) + 
        (Fsim42[(t-300),] %*% sd_fac[1:4,1:4] + t(mean_fac[1:4])) %*% t(Lsim42 %*% sd_lam[1:4,1:4] + rep(1,N) %*% t(mean_lam[1:4])) + 
        t(sqrt(diag(Ve[,4])) %*% rnorm(N,0,1))
    for (t in 401:500)
      Xsim[t,] = t(meanX) + 
        matrix(Fsim12[(t-400)]*sd_fac[1,1] + t(mean_fac[1]),1,1) %*% t(Lsim12*sd_lam[1,1] + rep(1,N)*mean_lam[1]) + 
        t(sqrt(diag(Ve[,1])) %*% rnorm(N,0,1))
    # stop simulation
    
    TT = as.matrix(c(100,120,200,220,300,320,400,420,500),9,1) # 9
    #                 1   4   4    8   8  4   4   1   1
    Tlocal = 25
    out = matrix(NA,10,ncol)
    
    kmax = 10
    eps1 = 1e-8
    eps2 = 0.25*c(0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80) # 0.68 good for Tshalf
    for (j in 1:9)
    {
   #   V = function(Xsim[(TT[j]-24):TT[j],],k) 
      temp2=matrix(NA,kmax,ncol)
      for (k in 1:kmax) 
        temp2[k,] = criteria(t(Xsim[(TT[j]-24):TT[j],]), k, kmax, eps1, eps2)
      l=1
      for (l in 1:ncol)
      {
        out[j,l] = which.min(temp2[,l]) 
        
      }
    }
    s=1
    for (s in 1:8)
    {    
      out[10,(12+4*s-3):(12+4*s)]  = rep(1,4)*eps2[s]
    }
  out
  outt = outt + out
  }  

  outt=outt/M
  }

####################################  
# righe 235-345 parallelizzate  1/4
simulations_parallel = function(mm, M, N, T,
                                meanX, mean_fac, mean_lam, sd_fac, sd_lam, Ve,
                                Fsim1, Fsim2, Fsim4, Fsim5, Fsim8, Fsim42, Fsim22, Fsim12,
                                Lsim1, Lsim2, Lsim4, Lsim5, Lsim8, Lsim22, Lsim42, Lsim12)
{
  set.seed(sample(1:M)[mm])
  ss1=sqrt(sd_fac[1,1]^2*sd_lam[1,1]^2/mean(Ve[,1]))
  ss4=sqrt(4*mean(diag(sd_fac[1:4,1:4]^2*sd_lam[1:4,1:4]^2))/mean(Ve[,4]))
  ss8=sqrt(8*mean(diag(sd_fac[1:8,1:8]^2*sd_lam[1:8,1:8]^2))/mean(Ve[,8]))
  
  # generazione dati
  Xsim = matrix(NA, nrow = T, ncol = N)
  for (t in 1:100)
    Xsim[t,] = t(meanX) + 
    (sd_fac[1,1]*Fsim1[t] + mean_fac[1]) %*% t(Lsim1*sd_lam[1,1] + mean_lam[1]) + 
    t(sqrt(diag(Ve[,1])) %*% rnorm(N,0,ss1))
  for (t in 101:200)
    Xsim[t,] = t(meanX) + 
    (Fsim4[(t-100),] %*% sd_fac[1:4,1:4] + t(mean_fac[1:4])) %*% t(Lsim4 %*% sd_lam[1:4,1:4] + rep(1,N) %*% t(mean_lam[1:4])) + 
    t(sqrt(diag(Ve[,4])) %*% rnorm(N,0,ss4))
  for (t in 201:300)
    Xsim[t,] = t(meanX) + 
    (Fsim8[(t-200),] %*% sd_fac[1:8,1:8] + t(mean_fac[1:8])) %*% t(Lsim8 %*% sd_lam[1:8,1:8] + rep(1,N) %*% t(mean_lam[1:8])) + 
    t(sqrt(diag(Ve[,8])) %*% rnorm(N,0,ss8))
  for (t in 301:400)
    Xsim[t,] = t(meanX) + 
    (Fsim42[(t-300),] %*% sd_fac[1:4,1:4] + t(mean_fac[1:4])) %*% t(Lsim42 %*% sd_lam[1:4,1:4] + rep(1,N) %*% t(mean_lam[1:4])) + 
    t(sqrt(diag(Ve[,4])) %*% rnorm(N,0,ss4))
  for (t in 401:500)
    Xsim[t,] = t(meanX) + 
    (sd_fac[1,1]*Fsim12[(t-400)] + mean_fac[1]) %*% t(Lsim12*sd_lam[1,1] + mean_lam[1]) + 
    t(sqrt(diag(Ve[,1])) %*% rnorm(N,0,ss1))
  # stop generazione dati
  
  TT = as.matrix(c(100,120,200,220,300,320,400,420,500),9,1) # 9
  #                 1   4   4    8   8  4   4   1   1
  Tlocal = 25
  
  # parametri stima
  eps1 = 1e-8
  eps2 = 0.25*c(0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80) # 0.68 good for Tshalf
  # eps2 = 0.65*c(0.25, 0.45) # 0.75 ok
  kmax = 10
  ncriteri = 12 + 4*length(eps2) + 1
  
  out  = matrix(NA, nrow = length(TT) + 1, ncol = ncriteri)
    for (j in 1:9)
    {
      # V = function(Xsim[(TT[j]-24):TT[j],],k) 
      temp2 = matrix(NA, nrow = kmax, ncol = ncriteri)
      for (k in 1:kmax) 
        temp2[k,] = criteria(t(Xsim[(TT[j]-24):TT[j],]), k, kmax, eps1, eps2)
      for (l in 1:ncriteri)
        out[j,l] = which.min(temp2[,l])
    }
  for (s in 1:8)
    out[10, (12 + 4*s - 3):(12 + 4*s)]  = rep(1,4)*eps2[s]
  return(out)
}
  

T = 500
N = 807#1243#807#283
set.seed(111208)
# genero Factors e Loadings
Fsim1  = matrix(rnorm(100,0,1), nrow = 100, ncol = 1)
Fsim2  = matrix(rnorm(200,0,1), nrow = 100, ncol = 2)
Fsim4  = matrix(rnorm(400,0,1), nrow = 100, ncol = 4)
Fsim5  = matrix(rnorm(500,0,1), nrow = 100, ncol = 5)
Fsim8  = matrix(rnorm(800,0,1), nrow = 100, ncol = 8)
Fsim42 = matrix(rnorm(400,0,1), nrow = 100, ncol = 4)
Fsim22 = matrix(rnorm(200,0,1), nrow = 100, ncol = 2)
Fsim12 = matrix(rnorm(100,0,1), nrow = 100, ncol = 1)

Lsim1  = matrix(rnorm(N*1,0,1), nrow = N, ncol = 1)
Lsim2  = matrix(rnorm(N*2,0,1), nrow = N, ncol = 2)
Lsim4  = matrix(rnorm(N*4,0,1), nrow = N, ncol = 4)
Lsim5  = matrix(rnorm(N*5,0,1), nrow = N, ncol = 5)
Lsim8  = matrix(rnorm(N*8,0,1), nrow = N, ncol = 8)
Lsim22 = matrix(rnorm(N*2,0,1), nrow = N, ncol = 2)
Lsim42 = matrix(rnorm(N*4,0,1), nrow = N, ncol = 4)
Lsim12 = matrix(rnorm(N*1,0,1), nrow = N, ncol = 1)
  
M = 100
supp = parLapply(cl = cl, X = 1:M, fun = simulations_parallel, 
                 M = M, N = N, T = T,
                 meanX = meanX, mean_fac = mean_fac, mean_lam = mean_lam,
                 sd_fac = sd_fac, sd_lam = sd_lam, Ve = Ve,
                 Fsim1 = Fsim1, Fsim2 = Fsim2, Fsim4 = Fsim4, Fsim5 = Fsim5, Fsim8 = Fsim8, Fsim42 = Fsim42, Fsim22 = Fsim22, Fsim12 = Fsim12,
                 Lsim1 = Lsim1, Lsim2 = Lsim2, Lsim4 = Lsim4, Lsim5 = Lsim5, Lsim8 = Lsim8, Lsim22 = Lsim22, Lsim42 = Lsim42, Lsim12 = Lsim12)
output = array(NA, dim = c(9+1, 12+4*8+1, M))
for (i in 1:M)
  output[,,i] = supp[[i]]
apply(output, 1:2, mean)
  
# Fine parallel
#############################################################
# Stesso esercizio di prima ma adesso loadings non time varying  2/4
simulations_parallel = function(mm, M, 
                                meanX, mean_fac, mean_lam, sd_fac, sd_lam, Ve,
                                Fsim, Lsim)
{
  set.seed(sample(1:M)[mm])
  N = nrow(Lsim) 
  T = nrow(Fsim) 
  
  
  ss1=sqrt(sd_fac[1,1]^2*sd_lam[1,1]^2/mean(Ve[,1]))
  ss4=sqrt(4*mean(diag(sd_fac[1:4,1:4]^2*sd_lam[1:4,1:4]^2))/mean(Ve[,4]))
  ss8=sqrt(8*mean(diag(sd_fac[1:8,1:8]^2*sd_lam[1:8,1:8]^2))/mean(Ve[,8]))
  
  nfact = c(rep(1, 100), rep(4, 100), rep(8, 100), rep(4, 100), rep(1, 100))
  nfacts = c(ss1*rep(1, 100), ss4*rep(1, 100), ss8*rep(1, 100), ss4*rep(1, 100), ss1*rep(1, 100))
  
  
  TT = matrix(c(100, 120, 200, 220, 300, 320, 400, 420, 500), nrow = 9, ncol = 1)
  #              1    4    4    8    8    4    4    1    1
  Tlocal = 25
  
  # generazione dati
  
  
  Xsim = matrix(NA, nrow = T, ncol = N)
  for (t in 1:T)
  {
    index = 1:nfact[t]
    Xsim[t,] = meanX + 
      (as.matrix(Lsim[,index]) %*% as.matrix(sd_lam[index,index]) + rep(1,N) %*% t(mean_lam[index])) %*% t(Fsim[t, index] %*% sd_fac[index,index] + t(mean_fac[index])) + 
      (sqrt(diag(Ve[,nfact[t]])) %*% rnorm(N, 0, nfacts[t]))
  }
  # stop generazione dati
  
  # parametri stima
  eps1 = 1e-8
  eps2 = 0.25*c(0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80) # 0.68 good for Tshalf
  # eps2 = 0.65*c(0.25, 0.45) # 0.75 ok
  kmax = 10
  ncriteri = 12 + 4*length(eps2) + 1
  
  out  = matrix(NA, nrow = length(TT) + 1, ncol = ncriteri)
  for (j in 1:9)
  {
    # V = function(Xsim[(TT[j]-24):TT[j],],k) 
    temp2 = matrix(NA, nrow = kmax, ncol = ncriteri)
    for (k in 1:kmax) 
      temp2[k,] = criteria(t(Xsim[(TT[j]-24):TT[j],]), k, kmax, eps1, eps2)
    for (l in 1:ncriteri)
      out[j,l] = which.min(temp2[,l])
  }
  for (s in 1:8)
    out[10, (12 + 4*s - 3):(12 + 4*s)]  = rep(1,4)*eps2[s]
  return(out)
}


T = 500
N = 807#1243#807#283
kmax = 10
# genero Factors e Loadings
set.seed(111208)
Fsim  = matrix(rnorm(T*kmax, 0, 1), nrow = T, ncol = kmax)
Lsim  = matrix(rnorm(N*kmax, 0, 1), nrow = N, ncol = kmax)

M = 100
supp = parLapply(cl = cl, X = 1:M, fun = simulations_parallel, 
                 M = M, meanX = meanX, mean_fac = mean_fac, mean_lam = mean_lam,
                 sd_fac = sd_fac, sd_lam = sd_lam, Ve = Ve,
                 Fsim  = Fsim, Lsim  = Lsim)
output = array(NA, dim = c(9+1, 12+4*8+1, M))
for (i in 1:M)
  output[,,i] = supp[[i]]
apply(output, 1:2, mean)


#############################################################
#  loadings  time varying BUT resampled  3/4



####################################  
# righe 235-345 parallelizzate
simulations_parallel = function(mm, M, N, T,
                                meanX, mean_fac, mean_lam, sd_fac, sd_lam, Ve,
                                Fsim1, Fsim2, Fsim4, Fsim5, Fsim8, Fsim42, Fsim22, Fsim12)#,
#                                Lsim1, Lsim2, Lsim4, Lsim5, Lsim8, Lsim22, Lsim42, Lsim12)
{
    set.seed(sample(1:M)[mm])
   # ratio=1#0.5
  #  ss1=ratio*sqrt(sd_fac[1,1]^2*sd_lam[1,1]^2/mean(Ve[,1]))
  #  ss4=ratio*sqrt(4*mean(diag(sd_fac[1:4,1:4]^2*sd_lam[1:4,1:4]^2))/mean(Ve[,4]))
  #  ss8=ratio*sqrt(8*mean(diag(sd_fac[1:8,1:8]^2*sd_lam[1:8,1:8]^2))/mean(Ve[,8]))
  # ss1=1
  # ss4=1
  # ss8=1
    
   #ss1=sqrt(mean(diag(var( Csim[1:100,]))))
   #ss4=sqrt(mean(diag(var( Csim[101:200,]))))
   #ss8=sqrt(mean(diag(var( Csim[201:300,]))))
   #ss42=sqrt(mean(diag(var( Csim[301:400,]))))
   #ss12=sqrt(mean(diag(var( Csim[401:500,]))))
   
  ratio=1#0.5
   ss1=0.03629963 
   ss4=0.04239013 
   ss8= 0.0455206
   ss42=0.04462886 
   ss12=0.03603078
   
   
    Lsim1  = matrix(rnorm(N*1,0,1), nrow = N, ncol = 1)
    Lsim2  = matrix(rnorm(N*2,0,1), nrow = N, ncol = 2)
    Lsim4  = matrix(rnorm(N*4,0,1), nrow = N, ncol = 4)
    Lsim5  = matrix(rnorm(N*5,0,1), nrow = N, ncol = 5)
    Lsim8  = matrix(rnorm(N*8,0,1), nrow = N, ncol = 8)
    Lsim22 = matrix(rnorm(N*2,0,1), nrow = N, ncol = 2)
    Lsim42 = matrix(rnorm(N*4,0,1), nrow = N, ncol = 4)
    Lsim12 = matrix(rnorm(N*1,0,1), nrow = N, ncol = 1)
    
    
    # generazione dati
    Xsim = matrix(NA, nrow = T, ncol = N)
    Csim = matrix(NA, nrow = T, ncol = N)
    esim = matrix(NA, nrow = T, ncol = N)
    nomean=1
    {
    
     for (t in 1:100)
  {    Xsim[t,] = nomean*t(meanX) + 
      (sd_fac[1,1]*Fsim1[t] + mean_fac[1]) %*% t(Lsim1*sd_lam[1,1] + mean_lam[1]) + 
    #  t(sqrt(diag(Ve[,1])) %*%
          ratio*rnorm(N,0,ss1)
    Csim[t,] = 
      (sd_fac[1,1]*Fsim1[t] + mean_fac[1]) %*% t(Lsim1*sd_lam[1,1] + mean_lam[1]) 
    
    esim[t,] =  rnorm(N,0,ss1)
    
     }
    A=eigen(sd_fac)
    AA=A$values
    AA[1:4]=AA[1]*t(rep(1,1,4))
    sd_fac2=A$vectors%*%diag(AA)%*%t(A$vectors)
    for (t in 101:200)
  {    Xsim[t,] = nomean*t(meanX) + 
      (Fsim4[(t-100),] %*% sd_fac2[1:4,1:4] + t(mean_fac[1:4])) %*% t(Lsim4 %*% sd_lam[1:4,1:4] + rep(1,N) %*% t(mean_lam[1:4])) + 
    #  t(sqrt(diag(Ve[,4])) %*% 
    ratio*rnorm(N,0,ss4)
  
    Csim[t,] =
      (Fsim4[(t-100),] %*% sd_fac2[1:4,1:4] + t(mean_fac[1:4])) %*% t(Lsim4 %*% sd_lam[1:4,1:4] + rep(1,N) %*% t(mean_lam[1:4]))
   
    esim[t,] =   rnorm(N,0,ss4)
    
    } 
    
    
    #  A=eigen(cov(Csim[201:300,]))
    #  B=cov(Csim[201:300,])
    # Csim2=Csim[201:300,]%*%solve(sqrtm(B))
    # AA=A$values
    # AA[2:4]=AA[1]*rep(1,3,1)
    # sAA=AA
    # sAA[1:30]=sqrt(AA[1:30])
    # Csim3=Csim2%*%(A$vectors%*%diag(sAA)%*%t(A$vectors))

        # 
    A=eigen(sd_fac)
    AA=A$values
    AA[1:8]=AA[1]*t(rep(1,1,8))
    sd_fac2=A$vectors%*%diag(AA)%*%t(A$vectors)
      for (t in 201:300)
      {
      Xsim[t,] = nomean*t(meanX) + 
      (Fsim8[(t-200),] %*% sd_fac2[1:8,1:8] + t(mean_fac[1:8])) %*% t(Lsim8 %*% sd_lam[1:8,1:8] + rep(1,N) %*% t(mean_lam[1:8]))+
        ratio*rnorm(N,0,ss8)   #  t(sqrt(diag(Ve[,8])) %*%

    
  
    Csim[t,] = 
      (Fsim8[(t-200),] %*% sd_fac2[1:8,1:8] + t(mean_fac[1:8])) %*% t(Lsim8 %*% sd_lam[1:8,1:8] + rep(1,N) %*% t(mean_lam[1:8])) 
     
    esim[t,] =   rnorm(N,0,ss8)
    
    
     }
      
    A=eigen(sd_fac)
    AA=A$values
    AA[1:4]=AA[1]*t(rep(1,1,4))
    sd_fac2=A$vectors%*%diag(AA)%*%t(A$vectors)
    for (t in 301:400)
    {
      Xsim[t,] = nomean*t(meanX) + 
      (Fsim42[(t-300),] %*% sd_fac2[1:4,1:4] + t(mean_fac[1:4])) %*% t(Lsim42 %*% sd_lam[1:4,1:4] + rep(1,N) %*% t(mean_lam[1:4])) + 
    #  t(sqrt(diag(Ve[,4])) %*%
        ratio*rnorm(N,0,ss42)
  
    Csim[t,] = 
      (Fsim42[(t-300),] %*% sd_fac2[1:4,1:4] + t(mean_fac[1:4])) %*% t(Lsim42 %*% sd_lam[1:4,1:4] + rep(1,N) %*% t(mean_lam[1:4])) 
    
    esim[t,] =  rnorm(N,0,ss42)
    
    
    } 
    
   
    for (t in 401:500)
    {
      Xsim[t,] = nomean*t(meanX) + 
      (sd_fac[1,1]*Fsim12[(t-400)] + mean_fac[1]) %*% t(Lsim12*sd_lam[1,1] + mean_lam[1]) + 
      #t(sqrt(diag(Ve[,1])) %*%
        ratio*rnorm(N,0,ss12)
   
    Csim[t,] = 
      (sd_fac[1,1]*Fsim12[(t-400)] + mean_fac[1]) %*% t(Lsim12*sd_lam[1,1] + mean_lam[1]) 
    
    esim[t,] =  rnorm(N,0,ss12)
    
    
    
    }
    }
    # # check
    # A=eigen(cov(Xsim[1:100,]))
    # A$values[1:10]
    # A=eigen(cov(Csim[1:100,]))
    # A$values[1:10]
    # 
    # 
    # 
    # A=eigen(cov(Xsim[101:200,]))
    # A$values[1:10]
    # A=eigen(cov(Csim[101:200,]))
    # A$values[1:10]
    # 
    # 
    # A=eigen(cov(Xsim[201:300,]))
    # A$values[1:15]
    # A=eigen(cov(Csim[201:300,]))
    # A$values[1:15]
    # 
    # 
    #     A=eigen(cov(Xsim[301:400,]))
    # A$values[1:15]
    # A=eigen(cov(Xsim[401:500,]))
    # A$values[1:15]
    
    
    # # check
    # 
   # mean(diag(var( Xsim[1:100,])))
  #  mean(diag(var( esim[1:100,])))
  #  mean(diag(var( Csim[1:100,])))
    # 
    # 
  #   mean(diag(var( Csim[1:100,])))/mean(diag(var( esim[1:100,])))
    # 
    # 
    #   mean(diag(var( Xsim[101:200,])))
    # mean(diag(var( Csim[101:200,])))
    # mean(diag(var( esim[101:200,])))
    # 
    # mean(diag(var( Csim[101:200,])))/mean(diag(var( esim[101:200,])))
    # 
    # 
    #     mean(diag(var( Xsim[201:300,])))
    # mean(diag(var( Csim[201:300,])))
    # mean(diag(var( esim[201:300,])))
    # 
    # mean(diag(var( Csim[201:300,])))/mean(diag(var( esim[201:300,])))
    # 
    #      
    # mean(diag(var( Xsim[301:400,])))
    # mean(diag(var( Csim[301:400,])))
    # mean(diag(var( esim[301:400,])))
    # 
    # mean(diag(var( Csim[301:400,])))/mean(diag(var( esim[301:400,])))
    # 
    #     
    # mean(diag(var( Xsim[401:500,])))
    # mean(diag(var( Csim[401:500,])))
    # mean(diag(var( esim[401:500,])))
    # 
    # mean(diag(var( Csim[401:500,])))/mean(diag(var( esim[401:500,])))
    # 
    # 
    # 
  
  
   # stop generazione dati
    Tlocal = 100 #25 #50#  10# 15  # divible by 4!
    
  #TT = as.matrix(c(100,(100+Tlocal*0.75),200,(200+Tlocal*0.75),300,(300+Tlocal*0.75),400,(400+Tlocal*0.75),500),9,1) # 9
  #                 1   4   4    8  8    8   4  4   1
    
    TT = as.matrix(c(100,200,300,400,500),5,1) # 9
  #                  1  4   8    4   1  
  
  # parametri stima
  eps1 = 1e-8
#eps2 = 0.25*c(0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80) # 0.68 good for Tshalf
 #  eps2 = 8*c(0.125, 0.135, 0.145, 0.155 , 0.165) # 0.75 ok
 #  eps2 = c(0.21,0.215, 0.22 ,  0.225, 0.23) # 0.75 ok
  # eps2 = c(0.235,0.24, 0.245 ,  0.25, 0.255) # 0.75 ok
 #eps2 = c(0.025,0.05, 0.075 ,  0.1, 0.125) # 0.75 ok
 # eps2 = c(-0.75,-0.6, -0.45 ,  -0.3, -0.15) # 0.75 ok
  
  # eps2 = c(-0.2,-0.15, -0.1 ,  -0.05, 0) # 0.75 ok
#PERFECT FOR MY CRITERIA AND T SMALL  eps2 = c(-0.8, -0.775 ,  -0.75, -0.725, -0.7) # 0.75 ok
 # eps2 = c(-0.3, -0.275 ,  -0.25, -0.225, -0.22) # perfect for T=10 (with sqrt(N) in penalty)
  eps2 = c(0.05, 0.1 ,  0.15 , 0.2 ,  0.25) # perfect for T=10 (with sqrt(N) in penalty)
  
  
  
  kmax = 10
  ncriteri = 12 + 4*length(eps2) + 1
  
  out  = matrix(NA, nrow = length(TT) + 2, ncol = ncriteri)
  for (j in 1:length(TT))
  {
    # V = function(Xsim[(TT[j]-24):TT[j],],k) 
    temp2 = matrix(NA, nrow = kmax, ncol = ncriteri)
    for (k in 1:kmax) 
      temp2[k,] = criteria(t(Xsim[(TT[j]-(Tlocal-1)):TT[j],]), k, kmax, eps1, eps2)
    for (l in 1:ncriteri)
      out[j,l] = which.min(temp2[,l])
  }
  for (s in 1:5)
    out[length(TT) + 1, (12 + 4*s - 3):(12 + 4*s)]  = rep(1,4)*eps2[s]
  
  out[length(TT) + 2,]=Tlocal*rep(1,1,33)
  out
  
  return(out)
}

TT = as.matrix(c(100,200,300,400,500),5,1) # 9
T = 500
N = 807#1243#807#283
set.seed(111208)
# genero Factors e Loadings
Fsim1  = matrix(rnorm(100,0,1), nrow = 100, ncol = 1)
Fsim2  = matrix(rnorm(200,0,1), nrow = 100, ncol = 2)
Fsim4  = matrix(rnorm(400,0,1), nrow = 100, ncol = 4)
Fsim5  = matrix(rnorm(500,0,1), nrow = 100, ncol = 5)
Fsim8  = matrix(rnorm(800,0,1), nrow = 100, ncol = 8)
Fsim42 = matrix(rnorm(400,0,1), nrow = 100, ncol = 4)
Fsim22 = matrix(rnorm(200,0,1), nrow = 100, ncol = 2)
Fsim12 = matrix(rnorm(100,0,1), nrow = 100, ncol = 1)


M = 50
supp = parLapply(cl = cl, X = 1:M, fun = simulations_parallel, 
                 M = M, N = N, T = T,
                 meanX = meanX, mean_fac = mean_fac, mean_lam = mean_lam,
                 sd_fac = sd_fac, sd_lam = sd_lam, Ve = Ve,
                 Fsim1 = Fsim1, Fsim2 = Fsim2, Fsim4 = Fsim4, Fsim5 = Fsim5, Fsim8 = Fsim8, Fsim42 = Fsim42, Fsim22 = Fsim22, Fsim12 = Fsim12)#,
#                 Lsim1 = Lsim1, Lsim2 = Lsim2, Lsim4 = Lsim4, Lsim5 = Lsim5, Lsim8 = Lsim8, Lsim22 = Lsim22, Lsim42 = Lsim42, Lsim12 = Lsim12)
output = array(NA, dim = c(length(TT) + 2, 12+5*4+1, M))
for (i in 1:M)
  output[,,i] = supp[[i]]
apply(output, 1:2, mean)



##################################################################



#############################################################
# Stesso esercizio di prima ma adesso loadings non time varying BUT resampled  4/4
simulations_parallel = function(mm, M, 
                                meanX, mean_fac, mean_lam, sd_fac, sd_lam, Ve,
                                Fsim,N)#, Lsim)
{
  set.seed(sample(1:M)[mm])
  
  Lsim  = matrix(rnorm(N*kmax, 0, 1), nrow = N, ncol = kmax)
  
  N = nrow(Lsim) 
  T = nrow(Fsim) 
  
  
  ss1=sqrt(sd_fac[1,1]^2*sd_lam[1,1]^2/mean(Ve[,1]))
  ss4=sqrt(4*mean(diag(sd_fac[1:4,1:4]^2*sd_lam[1:4,1:4]^2))/mean(Ve[,4]))
  ss8=sqrt(8*mean(diag(sd_fac[1:8,1:8]^2*sd_lam[1:8,1:8]^2))/mean(Ve[,8]))
  
  nfact = c(rep(1, 100), rep(4, 100), rep(8, 100), rep(4, 100), rep(1, 100))
  nfacts = c(ss1*rep(1, 100), ss4*rep(1, 100), ss8*rep(1, 100), ss4*rep(1, 100), ss1*rep(1, 100))
  
  
  TT = matrix(c(100, 120, 200, 220, 300, 320, 400, 420, 500), nrow = 9, ncol = 1)
  #              1    4    4    8    8    4    4    1    1
  Tlocal = 25
  
  # generazione dati
  
  
  Xsim = matrix(NA, nrow = T, ncol = N)
  for (t in 1:T)
  {
    index = 1:nfact[t]
    Xsim[t,] = meanX + 
      (as.matrix(Lsim[,index]) %*% as.matrix(sd_lam[index,index]) + rep(1,N) %*% t(mean_lam[index])) %*% t(Fsim[t, index] %*% sd_fac[index,index] + t(mean_fac[index])) + 
      (sqrt(diag(Ve[,nfact[t]])) %*% rnorm(N, 0, nfacts[t]))
  }
  # stop generazione dati
  
  # parametri stima
  eps1 = 1e-8
  eps2 = 0.25*c(0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80) # 0.68 good for Tshalf
  # eps2 = 0.65*c(0.25, 0.45) # 0.75 ok
  kmax = 10
  ncriteri = 12 + 4*length(eps2) + 1
  
  out  = matrix(NA, nrow = length(TT) + 1, ncol = ncriteri)
  for (j in 1:9)
  {
    # V = function(Xsim[(TT[j]-24):TT[j],],k) 
    temp2 = matrix(NA, nrow = kmax, ncol = ncriteri)
    for (k in 1:kmax) 
      temp2[k,] = criteria(t(Xsim[(TT[j]-24):TT[j],]), k, kmax, eps1, eps2)
    for (l in 1:ncriteri)
      out[j,l] = which.min(temp2[,l])
  }
  for (s in 1:8)
    out[10, (12 + 4*s - 3):(12 + 4*s)]  = rep(1,4)*eps2[s]
  return(out)
}


T = 500
N = 807#1243#807#283
kmax = 10
# genero Factors e Loadings
set.seed(111208)
Fsim  = matrix(rnorm(T*kmax, 0, 1), nrow = T, ncol = kmax)
#Lsim  = matrix(rnorm(N*kmax, 0, 1), nrow = N, ncol = kmax)

M = 100
supp = parLapply(cl = cl, X = 1:M, fun = simulations_parallel, 
                 M = M, meanX = meanX, mean_fac = mean_fac, mean_lam = mean_lam,
                 sd_fac = sd_fac, sd_lam = sd_lam, Ve = Ve,
                 Fsim  = Fsim, N=N)#, Lsim  = Lsim)
output = array(NA, dim = c(9+1, 12+4*8+1, M))
for (i in 1:M)
  output[,,i] = supp[[i]]
apply(output, 1:2, mean)

#  back to BaiNg DGP

####################################  
# righe 235-345 parallelizzate  5/4  risimulati & time varying
simulations_parallel = function(mm, M, N, T,
                                meanX, mean_fac, mean_lam, sd_fac, sd_lam, Ve,
                                Fsim1, Fsim2, Fsim4, Fsim5, Fsim8, Fsim42, Fsim22, Fsim12)#,
#                                Lsim1, Lsim2, Lsim4, Lsim5, Lsim8, Lsim22, Lsim42, Lsim12)
{
  set.seed(sample(1:M)[mm])
  ss1=1#sqrt(sd_fac[1,1]^2*sd_lam[1,1]^2/mean(Ve[,1]))
  ss4=sqrt(4)#sqrt(4*mean(diag(sd_fac[1:4,1:4]^2*sd_lam[1:4,1:4]^2))/mean(Ve[,4]))
  ss8=sqrt(8)#sqrt(8*mean(diag(sd_fac[1:8,1:8]^2*sd_lam[1:8,1:8]^2))/mean(Ve[,8]))
  
  Lsim1  = matrix(rnorm(N*1,0,1), nrow = N, ncol = 1)
  Lsim2  = matrix(rnorm(N*2,0,1), nrow = N, ncol = 2)
  Lsim4  = matrix(rnorm(N*4,0,1), nrow = N, ncol = 4)
  Lsim5  = matrix(rnorm(N*5,0,1), nrow = N, ncol = 5)
  Lsim8  = matrix(rnorm(N*8,0,1), nrow = N, ncol = 8)
  Lsim22 = matrix(rnorm(N*2,0,1), nrow = N, ncol = 2)
  Lsim42 = matrix(rnorm(N*4,0,1), nrow = N, ncol = 4)
  Lsim12 = matrix(rnorm(N*1,0,1), nrow = N, ncol = 1)
  
  ratio=0.1
  
  # generazione dati
  Xsim = matrix(NA, nrow = T, ncol = N)
  for (t in 1:100)
    Xsim[t,] = 
    Fsim1[t]%*% t(Lsim1) +  ratio*t(rnorm(N,0,ss1))
  for (t in 101:200)
    Xsim[t,] =  
    Fsim4[(t-100),]%*% t(Lsim4) + ratio*t(rnorm(N,0,ss4))
  for (t in 201:300)
    Xsim[t,] = 
    Fsim8[(t-200),] %*% t(Lsim8) + 
    ratio*t(rnorm(N,0,ss8))
  for (t in 301:400)
    Xsim[t,] =
    Fsim42[(t-300),]%*% t(Lsim42) + 
    ratio*t(rnorm(N,0,ss4))
  for (t in 401:500)
    Xsim[t,] = 
Fsim12[(t-400)]%*% t(Lsim12) + 
    ratio*t(rnorm(N,0,ss1))
  # stop generazione dati
  
  # # check
  # A=eigen(cov(Xsim[1:100,]))
  # A$values[1:10]
  # A=eigen(cov(Csim[1:100,]))
  # A$values[1:10]
  # 
  # 
  # 
  # A=eigen(cov(Xsim[101:200,]))
  # A$values[1:10]
  # A=eigen(cov(Csim[101:200,]))
  # A$values[1:10]
  # 
  # 
  # A=eigen(cov(Xsim[201:300,]))
  # A$values[1:15]
  # A=eigen(cov(Csim[201:300,]))
  # A$values[1:15]
  # 
  # 
  #     A=eigen(cov(Xsim[301:400,]))
  # A$values[1:15]
  # A=eigen(cov(Xsim[401:500,]))
  # A$values[1:15]
  
  
  # # check
  # 
  
  
  
  TT = as.matrix(c(100,200,300,400,500),5,1) # 9
  #                 1   4   8    4   1  
  Tlocal = 15
  
  # parametri stima
  eps1 = 1e-8
 # eps2 = 0.25*c(0.4, 1, 1.8, 0.60, 0.65, 0.70, 0.75, 0.80) # 0.68 good for Tshalf
 eps2 = c(0.1, 0.25, 0.45) # 0.68 good for Tshalf
  
    kmax = 10
  ncriteri = 12 + 4*length(eps2) + 1
  
  out  = matrix(NA, nrow = length(TT) + 2, ncol = ncriteri)
  for (j in 1:5)
  {
    # V = function(Xsim[(TT[j]-24):TT[j],],k) 
    temp2 = matrix(NA, nrow = kmax, ncol = ncriteri)
    for (k in 1:kmax) 
      temp2[k,] = criteria(t(Xsim[(TT[j]-(Tlocal-1)):TT[j],]), k, kmax, eps1, eps2)
    for (l in 1:ncriteri)
      out[j,l] = which.min(temp2[,l])
  }
  for (s in 1:3)
    out[6, (12 + 4*s - 3):(12 + 4*s)]  = rep(1,4)*eps2[s]
  out[7,]=Tlocal*rep(1,1,25) 
  
    return(out)
}


T = 500
N = 807#1243#807#283
set.seed(111208)
# genero Factors e Loadings
Fsim1  = matrix(rnorm(100,0,1), nrow = 100, ncol = 1)
Fsim2  = matrix(rnorm(200,0,1), nrow = 100, ncol = 2)
Fsim4  = matrix(rnorm(400,0,1), nrow = 100, ncol = 4)
Fsim5  = matrix(rnorm(500,0,1), nrow = 100, ncol = 5)
Fsim8  = matrix(rnorm(800,0,1), nrow = 100, ncol = 8)
Fsim42 = matrix(rnorm(400,0,1), nrow = 100, ncol = 4)
Fsim22 = matrix(rnorm(200,0,1), nrow = 100, ncol = 2)
Fsim12 = matrix(rnorm(100,0,1), nrow = 100, ncol = 1)

# Lsim1  = matrix(rnorm(N*1,0,1), nrow = N, ncol = 1)
# Lsim2  = matrix(rnorm(N*2,0,1), nrow = N, ncol = 2)
# Lsim4  = matrix(rnorm(N*4,0,1), nrow = N, ncol = 4)
# Lsim5  = matrix(rnorm(N*5,0,1), nrow = N, ncol = 5)
# Lsim8  = matrix(rnorm(N*8,0,1), nrow = N, ncol = 8)
# Lsim22 = matrix(rnorm(N*2,0,1), nrow = N, ncol = 2)
# Lsim42 = matrix(rnorm(N*4,0,1), nrow = N, ncol = 4)
# Lsim12 = matrix(rnorm(N*1,0,1), nrow = N, ncol = 1)

M = 100
supp = parLapply(cl = cl, X = 1:M, fun = simulations_parallel, 
                 M = M, N = N, T = T,
                 meanX = meanX, mean_fac = mean_fac, mean_lam = mean_lam,
                  sd_fac = sd_fac, sd_lam = sd_lam, Ve = Ve,
                 Fsim1 = Fsim1, Fsim2 = Fsim2, Fsim4 = Fsim4, Fsim5 = Fsim5, Fsim8 = Fsim8, Fsim42 = Fsim42, Fsim22 = Fsim22, Fsim12 = Fsim12)#,
#                 Lsim1 = Lsim1, Lsim2 = Lsim2, Lsim4 = Lsim4, Lsim5 = Lsim5, Lsim8 = Lsim8, Lsim22 = Lsim22, Lsim42 = Lsim42, Lsim12 = Lsim12)
output = array(NA, dim = c(6+1, 12+4*3+1, M))
for (i in 1:M)
  output[,,i] = supp[[i]]
apply(output, 1:2, mean)

# Fine parallel
#############################################################
# Stesso esercizio di prima ma adesso loadings non time varying  6/4
simulations_parallel = function(mm, M, 
                                meanX, mean_fac, mean_lam, sd_fac, sd_lam, Ve,
                                Fsim, N,kmax)#, Lsim)
{
  set.seed(sample(1:M)[mm])
  #N = nrow(Lsim) 
  T = nrow(Fsim) 
  
  Lsim  = matrix(rnorm(N*kmax, 0, 1), nrow = N, ncol = kmax)
  
  ss1=1#sqrt(sd_fac[1,1]^2*sd_lam[1,1]^2/mean(Ve[,1]))
  ss4=2#sqrt(4*mean(diag(sd_fac[1:4,1:4]^2*sd_lam[1:4,1:4]^2))/mean(Ve[,4]))
  ss8=sqrt(8)#sqrt(8*mean(diag(sd_fac[1:8,1:8]^2*sd_lam[1:8,1:8]^2))/mean(Ve[,8]))
  
  nfact = c(rep(1, 100), rep(4, 100), rep(8, 100), rep(4, 100), rep(1, 100))
  nfacts = c(ss1*rep(1, 100), ss4*rep(1, 100), ss8*rep(1, 100), ss4*rep(1, 100), ss1*rep(1, 100))
  
  
  TT = matrix(c(100, 120, 200, 220, 300, 320, 400, 420, 500), nrow = 9, ncol = 1)
  #              1    4    4    8    8    4    4    1    1
  Tlocal = 25
  
  # generazione dati
  
  
  Xsim = matrix(NA, nrow = T, ncol = N)
  for (t in 1:T)
  {
    index = 1:nfact[t]
    Xsim[t,] = meanX + 
      (as.matrix(Lsim[,index]) %*% as.matrix(sd_lam[index,index]) + rep(1,N) %*% t(mean_lam[index])) %*% t(Fsim[t, index] %*% sd_fac[index,index] + t(mean_fac[index])) + 
      (sqrt(diag(Ve[,nfact[t]])) %*% rnorm(N, 0, nfacts[t]))
  }
  # stop generazione dati
  
  # parametri stima
  eps1 = 1e-8
  #eps2 = 0.25*c(0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80) # 0.68 good for Tshalf
   eps2 = c(0.1 ,0.25, 0.45) # 0.75 ok
  kmax = 10
  ncriteri = 12 + 4*length(eps2) + 1
  
  out  = matrix(NA, nrow = length(TT) + 1, ncol = ncriteri)
  for (j in 1:9)
  {
    # V = function(Xsim[(TT[j]-24):TT[j],],k) 
    temp2 = matrix(NA, nrow = kmax, ncol = ncriteri)
    for (k in 1:kmax) 
      temp2[k,] = criteria(t(Xsim[(TT[j]-24):TT[j],]), k, kmax, eps1, eps2)
    for (l in 1:ncriteri)
      out[j,l] = which.min(temp2[,l])
  }
  for (s in 1:3)
    out[10, (12 + 4*s - 3):(12 + 4*s)]  = rep(1,4)*eps2[s]
  return(out)
}


T = 500
N = 807#1243#807#283
kmax = 10
# genero Factors e Loadings
set.seed(111208)
Fsim  = matrix(rnorm(T*kmax, 0, 1), nrow = T, ncol = kmax)
#Lsim  = matrix(rnorm(N*kmax, 0, 1), nrow = N, ncol = kmax)

M = 100
supp = parLapply(cl = cl, X = 1:M, fun = simulations_parallel, 
                 M = M, meanX = meanX, mean_fac = mean_fac, mean_lam = mean_lam,
                 sd_fac = sd_fac, sd_lam = sd_lam, Ve = Ve,
                 Fsim  = Fsim, N=N, kmax=kmax)# Lsim  = Lsim)
output = array(NA, dim = c(9+1, 12+4*3+1, M))
for (i in 1:M)
  output[,,i] = supp[[i]]
apply(output, 1:2, mean)


#############################################################
#


#STOP  TIME VARYING  NUMBER OF FACTORS

# START CONSTANT NUMBER OF FACTORS

#  heterosk.ts = NA
#  corr.ts = NA
#  corr.cs = NA
  
  
  
  # mean_lam = apply(lambda_dgp,2,mean)
  # var_lam = cov(lambda_dgp)
  # 
  # mean_fac = apply(factors_dgp,2,mean)
  # var_fac = cov(factors_dgp)
  # 
  # e1=Xd-factors_dgp[,1]%*%t(lambda_dgp[,1])
  # V1=diag(diag(cov(e1)))
  
  
  set.seed(12)
  table2 = tabella2(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table2, file = "Table_I_add.csv")
  # extensions
  N.vec = c(3000, 9000)
  T.vec = c(5, 5)
  
  set.seed(12)
  table2ext = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table2ext, file = "Table_I_add_ext.csv")
  
  # T very small 
  kmax = 2
  n.iter = 250
  N.vec = c(5, 15, 30, 100, 200, 500, 1000, 3000)
  T.vec = rep(2, length(N.vec))
  set.seed(12)
  table1 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table1, file = "1fac_T=2")
  
  kmax = 3
  n.iter = 250
  N.vec = c(5, 15, 30, 100, 200, 500, 1000, 3000)
  T.vec = rep(3, length(N.vec))
  set.seed(12)
  table2 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table2, file = "1fac_T=3")
}





# Homoskedstic, uncorrelated, var(common signal) = var(error) (table I, II, III)
# 1 FATTORE (table I)
{
  # Replica tabella 1 BAI NG
  {
  theta = r = 3
  kmax = 8
  n.iter = 30#250
  N.vec = c(100)
  T.vec = c(5)
  
  #N.vec = c(100, 100, 200, 500, 1000, 2000)
  #T.vec = c(40, 60, 60, 60, 60, 60)
  eps1 = NA
  eps2 = NA
  heterosk.ts = NA
  corr.ts = NA
  corr.cs = NA
  
  set.seed(12)
  table1 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table1, file = "Table_I.csv")
  }
  
  # tabella con dati richiesti da Prof
  {
  theta = r = 1
  kmax = 8
  n.iter = 7
#  N.vec = rep(c(10, 30, 100, 200, 500, 1000, 3000), 5)
#  T.vec = c(rep(5, length(unique(N.vec))), rep(10, length(unique(N.vec))), rep(30, length(unique(N.vec))), rep(60, length(unique(N.vec))), rep(120, length(unique(N.vec))))
    N.vec = c(10,20)#rep(c(10), 2)
    T.vec = c(rep(5, length(unique(N.vec))))
  
    eps1 = 1e-8
  eps2 = c(0.1, 0.25, 0.45)
  heterosk.ts = NA
  corr.ts = NA
  corr.cs = NA
  
  # mean_lam = apply(lambda_dgp,2,mean)
  # var_lam = cov(lambda_dgp)
  # 
  # mean_fac = apply(factors_dgp,2,mean)
  # var_fac = cov(factors_dgp)
  # 
  # e1=Xd-factors_dgp[,1]%*%t(lambda_dgp[,1])
  # V1=diag(diag(cov(e1)))
  
  
  set.seed(12)
  table2 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table2, file = "Table_I_add.csv")
  # extensions
  N.vec = c(3000, 9000)
  T.vec = c(5, 5)
  
  set.seed(12)
  table2ext = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table2ext, file = "Table_I_add_ext.csv")
  
  # T very small 
  kmax = 2
  n.iter = 250
  N.vec = c(5, 15, 30, 100, 200, 500, 1000, 3000)
  T.vec = rep(2, length(N.vec))
  set.seed(12)
  table1 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table1, file = "1fac_T=2")
  
  kmax = 3
  n.iter = 250
  N.vec = c(5, 15, 30, 100, 200, 500, 1000, 3000)
  T.vec = rep(3, length(N.vec))
  set.seed(12)
  table2 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table2, file = "1fac_T=3")
  }
  
  
  
}
# 3 FATTORI (table II)
{
  # Replica tabella 2 BAI NG
  {
  theta = r = 3
  kmax = 8
  n.iter = 250
  N.vec = c(100, 100, 200, 500, 1000, 2000)
  T.vec = c(40, 60, 60, 60, 60, 60)
  eps1 = NA
  eps2 = NA
  heterosk.ts = NA
  corr.ts = NA
  corr.cs = NA
  
  set.seed(12)
  table3 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table3, file = "Table_II.csv")
  }
  
  # tabella con dati richiesti da Prof
  {
  theta = r = 3
  kmax = 8
  n.iter = 100
  N.vec = rep(c(10, 30, 100, 200, 500, 1000, 3000), 5)
  T.vec = c(rep(5, length(unique(N.vec))), rep(10, length(unique(N.vec))), rep(30, length(unique(N.vec))), rep(60, length(unique(N.vec))), rep(120, length(unique(N.vec))))
  eps1 = 1e-8
  eps2 = c(0.1, 0.25, 0.45)
  heterosk.ts = NA
  corr.ts = NA
  corr.cs = NA
  
  set.seed(12)
  table4 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table4, file = "Table_II_add.csv")
  # extensions
  N.vec = c(3000, 9000)
  T.vec = c(5, 5)
  
  set.seed(12)
  table4ext = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table4ext, file = "Table_II_add_ext.csv")
  }
}
# 5 FATTORI (table III)
{
  # Replica tabella 3 BAI NG
  {
  theta = r = 5
  kmax = 8
  n.iter = 250
  N.vec = c(100, 100, 200, 500, 1000, 2000)
  T.vec = c(40, 60, 60, 60, 60, 60)
  eps1 = NA
  eps2 = NA
  heterosk.ts = NA
  corr.ts = NA
  corr.cs = NA
  
  set.seed(12)
  table5 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table5, file = "Table_III.csv")
  }
  
  # tabella con dati richiesti da Prof
  {
  theta = r = 5
  kmax = 8
  n.iter = 100
  N.vec = rep(c(10, 30, 100, 200, 500, 1000, 3000), 5)
  T.vec = c(rep(5, length(unique(N.vec))), rep(10, length(unique(N.vec))), rep(30, length(unique(N.vec))), rep(60, length(unique(N.vec))), rep(120, length(unique(N.vec))))
  eps1 = 1e-8
  eps2 = c(0.1, 0.25, 0.45)
  heterosk.ts = NA
  corr.ts = NA
  corr.cs = NA
  
  set.seed(12)
  table6 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table6, file = "Table_III_add.csv")
  # extensions
  N.vec = c(3000, 9000)
  T.vec = c(5, 5)
  
  set.seed(12)
  table6ext = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table6ext, file = "Table_III_add_ext.csv")
  
  
  N.vec = c(10, 30, 100, 200)#, 500, 1000, 3000, 5000)
  T.vec = rep(10, length(unique(N.vec)))
  
  kmax = 9
  set.seed(12)
  table_a = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  
  kmax = 10
  set.seed(12)
  table_aa = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  
  kmax = 15
  set.seed(12)
  table_aaa = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  }
}  
# 10 FATTORI
{
  # tabella con dati richiesti da Prof
  {
  theta = r = 10
  kmax = 15
  n.iter = 100
  N.vec = rep(c(15, 30, 100, 200, 500, 1000, 3000), 5)
  T.vec = c(rep(5, length(unique(N.vec))), rep(10, length(unique(N.vec))), rep(30, length(unique(N.vec))), rep(60, length(unique(N.vec))), rep(120, length(unique(N.vec))))
  eps1 = 1e-8
  eps2 = c(0.1, 0.25, 0.45)
  heterosk.ts = NA
  corr.ts = NA
  corr.cs = NA
  
  set.seed(12)
  table7 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table7, file = "Table_10fac.csv")
  
  # extensions
  N.vec = c(3000, 9000)
  T.vec = c(5, 5)
  
  set.seed(12)
  table7ext = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table7ext, file = "Table_10fac_ext.csv")
  }
}

# Homoskedstic, uncorrelated, var(common signal) < var(error) (table V)
# 5 FATTORI (table V)
{
  # Replica tabella 5 BAI NG
  r = 5
  theta = r*2
  kmax = 8
  n.iter = 250
  N.vec = c(100, 100, 200, 500, 1000, 2000)
  T.vec = c(40, 60, 60, 60, 60, 60)
  eps1 = NA
  eps2 = NA
  heterosk.ts = NA
  corr.ts = NA
  corr.cs = NA
  
  set.seed(12)
  table8 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table8, file = "Table_V.csv")
  
  # tabella con dati richiesti da Prof
  r = 5
  theta = r*2  
  kmax = 8
  n.iter = 100
  N.vec = rep(c(10, 30, 100, 200, 500, 1000, 3000), 5)
  T.vec = c(rep(5, length(unique(N.vec))), rep(10, length(unique(N.vec))), rep(30, length(unique(N.vec))), rep(60, length(unique(N.vec))), rep(120, length(unique(N.vec))))
  eps1 = 1e-8
  eps2 = c(0.1, 0.25, 0.45)
  heterosk.ts = NA
  corr.ts = NA
  corr.cs = NA
  
  set.seed(12)
  table9 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table9, file = "Table_V_add.csv")
}

# Heteroskedastic and uncorrelated (table IV)
# 5 FATTORI (table IV)
{
  # Replica tabella 3 BAI NG
  r = 5
  theta = r
  kmax = 8
  n.iter = 250
  N.vec = c(100, 100, 200, 500, 1000, 2000)
  T.vec = c(40, 60, 60, 60, 60, 60)
  eps1 = NA
  eps2 = NA
  heterosk.ts = T
  corr.ts = NA
  corr.cs = NA
  
  set.seed(12)
  table10 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table10, file = "Table_IV.csv")
  
  # tabella con dati richiesti da Prof
  r = 5
  theta = r
  kmax = 8
  n.iter = 100
  N.vec = rep(c(10, 30, 100, 200, 500, 1000, 3000), 5)
  T.vec = c(rep(5, length(unique(N.vec))), rep(10, length(unique(N.vec))), rep(30, length(unique(N.vec))), rep(60, length(unique(N.vec))), rep(120, length(unique(N.vec))))
  eps1 = 1e-8
  eps2 = c(0.1, 0.25, 0.45)
  heterosk.ts = T
  corr.ts = NA
  corr.cs = NA
  
  set.seed(12)
  table11 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table11, file = "Table_IV_add.csv")
}

# Serial and cross correlation (tables VI, VII, VIII)
# Only serial
# 5 FATTORI (table VI)
{
  # Replica tabella 6 BAI NG
  r = 5
  theta = r
  kmax = 8
  n.iter = 250
  N.vec = c(100, 100, 200, 500, 1000, 2000)
  T.vec = c(40, 60, 60, 60, 60, 60)
  eps1 = NA
  eps2 = NA
  heterosk.ts = NA
  corr.ts = 0.5
  corr.cs = NA
  
  set.seed(12)
  table12 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table12, file = "Table_VI.csv")
  
  # tabella con dati richiesti da Prof
  r = 5
  theta = r
  kmax = 8
  n.iter = 100
  N.vec = rep(c(10, 30, 100, 200, 500, 1000, 3000), 5)
  T.vec = c(rep(5, length(unique(N.vec))), rep(10, length(unique(N.vec))), rep(30, length(unique(N.vec))), rep(60, length(unique(N.vec))), rep(120, length(unique(N.vec))))
  eps1 = 1e-8
  eps2 = c(0.1, 0.25, 0.45)
  heterosk.ts = NA
  corr.ts = 0.5
  corr.cs = NA
  
  set.seed(12)
  table13 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table13, file = "Table_VI_add.csv")
}
# Only cross-section
# 5 FATTORI (table VII)
{
  # Replica tabella 7 BAI NG
  r = 5
  theta = r
  kmax = 8
  n.iter = 250
  N.vec = c(100, 100, 200, 500, 1000, 2000)
  T.vec = c(40, 60, 60, 60, 60, 60)
  eps1 = NA
  eps2 = NA
  heterosk.ts = NA
  corr.ts = NA
  corr.cs = 0.2
  
  set.seed(12)
  table14 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table14, file = "Table_VII.csv")
  
  # tabella con dati richiesti da Prof
  r = 5
  theta = r
  kmax = 8
  n.iter = 100
  N.vec = rep(c(10, 30, 100, 200, 500, 1000, 3000), 5)
  T.vec = c(rep(5, length(unique(N.vec))), rep(10, length(unique(N.vec))), rep(30, length(unique(N.vec))), rep(60, length(unique(N.vec))), rep(120, length(unique(N.vec))))
  eps1 = 1e-8
  eps2 = c(0.1, 0.25, 0.45)
  heterosk.ts = NA
  corr.ts = NA
  corr.cs = 0.2
  
  set.seed(12)
  table15 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table15, file = "Table_VII_add.csv")
}
# Both
# 5 FATTORI (table VIII)
{
  # Replica tabella 8 BAI NG
  r = 5
  theta = r
  kmax = 8
  n.iter = 250
  N.vec = c(100, 100, 200, 500, 1000, 2000)
  T.vec = c(40, 60, 60, 60, 60, 60)
  eps1 = NA
  eps2 = NA
  heterosk.ts = NA
  corr.ts = 0.5
  corr.cs = 0.2
  
  set.seed(12)
  table16 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table16, file = "Table_VIII.csv")
  
  # tabella con dati richiesti da Prof
  r = 5
  theta = r
  kmax = 8
  n.iter = 100
  N.vec = rep(c(10, 30, 100, 200, 500, 1000, 3000), 5)
  T.vec = c(rep(5, length(unique(N.vec))), rep(10, length(unique(N.vec))), rep(30, length(unique(N.vec))), rep(60, length(unique(N.vec))), rep(120, length(unique(N.vec))))
  eps1 = 1e-8
  eps2 = c(0.1, 0.25, 0.45)
  heterosk.ts = NA
  corr.ts = 0.5
  corr.cs = 0.2
  
  set.seed(12)
  table17 = tabella(r, theta, kmax, N.vec, T.vec, n.iter, eps1, eps2, heterosk.ts, corr.ts, corr.cs)
  write.csv(table17, file = "Table_VIII_add.csv")
}